/*This program calculates the Pearson's correlation coefficient as well as two
  rank correlation coefficient, namely, spearman's and kendall's

  author: Chao Wang
  Date: Sep. 2nd, 2012
    
  Input pairs are limited to the number of 1024
*/
          
#include<stdio.h>
#include<stdlib.h>
#include<assert.h>
#include<math.h>
#include<memory.h>

#define N 1024

//calculate the mean of the elements in array x of size n 
double mean(double x[], int n)
{
       double sum = 0.0f;
       int i = 0;
       for(i = 0; i < n; ++i)
             sum += x[i];
       return sum/n;       
}

//calculate the Pearson's correlation coefficient 
double calculate_pearsons_p(double x[], double y[], int n)
{
       double x_mean = mean(x, n);
       double y_mean = mean(y, n);
       
       double up = 0.0f, down1 = 0.0f, down2 = 0.0f;
       int i = 0;
       for(i = 0; i < n; ++i)
       {
             up += (x[i] - x_mean) * (y[i] - y_mean);
             down1 += (x[i] - x_mean) * (x[i] - x_mean);
             down2 += (y[i] - y_mean) * (y[i] - y_mean);
       }     
       
       return up/sqrt(down1*down2);
}


int bubblesort(double a[], int start, int n)
{
	int i, j;
	for(i = start; i < start + n; i++)
	{
		for(j = i + 1; j < start + n; j++)
		{
			if(a[i] > a[j])
			{
				double tmp;
				tmp = a[i];
				a[i] = a[j];
				a[j] = tmp;
			}
		}
	}
}


//binary searching an element in the array of size n
int binary_search(const double array[], const int n, const double element)
{
       int low = 1, high = n, mid;
       while(low <= high)
       {
            mid = (low + high)>>1;
            if(array[mid] - element > 1e-5)
                high = mid - 1;
            else if(array[mid] - element < -1e-5)
                low = mid + 1;
	    else
                return mid;
       }
       
       return -1;    
}

//calcuate the spearman's rank correlation coefficient
double calculate_spearman_rank_p(double x[], double y[], int n)
{
       //we only use the first n positions for holding the ranks of elements in array x and y.
       double rank_x[N] = {0};
       double rank_y[N] = {0};
       
       
       double * x_tmp = (double *)malloc(sizeof(double)*(n+1));
       double * y_tmp = (double *)malloc(sizeof(double)*(n+1));
       
       memcpy(x_tmp + 1, x, sizeof(double)*n);
       memcpy(y_tmp + 1, y, sizeof(double)*n);
       
       bubblesort(x_tmp, 1, n);
       bubblesort(y_tmp, 1, n);
       //using binary searching method to find the 'raw' ranks of each element
       int i = 0;
       for(i = 0; i < n; ++i)
       {
             int raw_rank = binary_search(x_tmp, n, x[i]);
             if(raw_rank != -1)
             {
                 //deal with tried elements in the array
                 int ll, hh, j;
                 for(ll = raw_rank; ll >= 1 && x_tmp[ll] == x_tmp[raw_rank]; ll--)
                         ;
                 for(hh = raw_rank; hh <= n && x_tmp[hh] == x_tmp[raw_rank]; hh++)
                         ;
                 
                 int ties_count = hh-ll-1, avg = 0;
                 ll++;
                 hh--;        
                 for(j = ll; j <= hh; j++)
                       avg += j;
                         
                 rank_x[i] = ((double)avg)/ties_count;
             }
             else
                 assert(0);      
       }

       for(i = 0; i < n; ++i)
       {
             int raw_rank = binary_search(y_tmp, n, y[i]);
             if(raw_rank != -1)
             {
                 //deal with tried elements in the array
                 int ll, hh, j;
                 for(ll = raw_rank; ll >= 1 && y_tmp[ll] == y_tmp[raw_rank]; ll--)
                         ;
                 for(hh = raw_rank; hh <= n && y_tmp[hh] == y_tmp[raw_rank]; hh++)
                         ;
                 
                 int ties_count = hh-ll-1, avg = 0;        
                 ll++;
                 hh--;
                 for(j = ll; j <= hh; j++)
                       avg += j;
                         
                 rank_y[i] = ((double)avg)/ties_count;
             }
             else
                 assert(0);      
       }
       
       double rank_x_mean = mean(rank_x, n);
       double rank_y_mean = mean(rank_y, n);
       
       double up = 0.0f, down1 = 0.0f, down2 = 0.0f;
       for(i = 0; i < n; ++i)
       {
             up += (rank_x[i] - rank_x_mean) * (rank_y[i] - rank_y_mean);
             down1 += (rank_x[i] - rank_x_mean) * (rank_x[i] - rank_x_mean);
             down2 += (rank_y[i] - rank_y_mean) * (rank_y[i] - rank_y_mean);
       }     
       
       free(x_tmp);
       free(y_tmp); 
       return up/sqrt(down1*down2);
}

//determine if two pairs (x1, y1) and (x2, y2) are concordant
int concordant(double x1, double x2, double y1, double y2)
{
    if(((x1 > x2) && (y1 > y2)) || ((x1 < x2) && (y1 < y2)))
            return 1;
    else
            return 0;                  
}

//determine if two pairs (x1, y1) and (x2, y2) are disconcordant
int disconcordant(double x1, double x2, double y1, double y2)
{
    if(((x1 > x2) && (y1 < y2)) || ((x1 < x2) && (y1 > y2)))
            return 1;
    else
            return 0;
}

//calculate the kendall rank correlation coefficient
double calculate_kendall_rank_p(double x[], double y[], int n)
{
       double n0 = n*(n-1)/2;
       double n1 = 0.0f;
       double n2 = 0.0f;
       double nc = 0.0f;
       double nd = 0.0f;
       
       double * x_tmp = (double*)malloc(sizeof(double)*n);
       double * y_tmp = (double*)malloc(sizeof(double)*n);
       
       memcpy(x_tmp, x, sizeof(double)*n);
       memcpy(y_tmp, y, sizeof(double)*n);
      
       bubblesort(x_tmp, 0, n);
       bubblesort(y_tmp, 0, n);
       
       int i = 0, j = 0, t = 0;
       for(i = 0; i < n; ++i)
       {
           for(j = i; j < n && x_tmp[j] == x_tmp[i]; ++j)
                 ;
           t = j - i;
           n1 += t*(t-1)/2;
           i = j - 1;
       }
       for(i = 0; i < n; ++i)
       {
           for(j = i; j < n && y_tmp[j] == y_tmp[i]; ++j)
                 ;
           t = j - i;
           n2 += t*(t-1)/2;
           i = j - 1;      
       }
              
       for(i = 0; i < n; ++i)
       {
             for(j = i + 1; j < n; ++j)
             {
                  if(concordant(x[i], x[j], y[i], y[j]))
                        nc++;
                  else if(disconcordant(x[i], x[j], y[i], y[j]))
                        nd++; 
             }
       }
       
       free(x_tmp);
       free(y_tmp);
       return (nc-nd)/sqrt((n0-n1)*(n0-n2));
}

int main()
{
    double x[N] = {0};
    double y[N] = {0};
    int n;
    
    printf("insert the number of paird trials:\n");
    scanf("%d", &n);
    
    int i = 0;
    for(i = 0; i < n; ++i)
    {
          scanf("%lf%lf", &x[i], &y[i]);
    }

    printf("pearsons' correlation coefficient is: %lf\n", calculate_pearsons_p(x,y,n));
    printf("spearman' correlation coefficient is: %lf\n", calculate_spearman_rank_p(x,y,n));
    printf("kendall' correlation coefficient is: %lf\n", calculate_kendall_rank_p(x,y,n));  
    
    return 0;
}
